  # Clearing environtment
  rm(list = ls())
  
  # Loading Packages
  library(tidyverse)
  library(magrittr)
  
  # Loading Data
  dat <- read.csv(file = "./RILA_Data_final.csv")
  
  # Save options
  sa <- 1
  width <- 32
  height <- 16
  
  width2 <- 36
  t1 <- 45
  t2 <- 45
  
  # Overleaf Path
  opath <- "../../../Apps/Overleaf/Target Date RILAs/"
  
  # Cleaning Data -----------------------------------------------------------
  # Removing too short
  dat %<>% 
    filter((finish - start)/60 > 30) %>% 
    filter(grepl("385", la10) &
             grepl("200", la11))
  
  # Retirement account only
  dat <- dat %>%  filter(aq6 == "Yes")
  
  # l <- which(table(dat$Worker_ID)>1)
  # k <- which(dat$Worker_ID %in% names(table(dat$Worker_ID))[l])
  # 
  # dat <- dat[-k,]
  
  cor(dat$group_assigned, dat$gross_expense_ratio_random)
  dat$change1 <- ifelse(dat$q1_a == dat$q1_b, 0, 1) 
  dat$change2 <- ifelse(dat$q2_a == dat$q2_b, 0, 1) 
  dat$change3 <- ifelse(dat$q3_a == dat$q3_b, 0, 1) 
  dat$change4 <- ifelse(dat$q4_a == dat$q4_b, 0, 1) 
  
  # Labeling Answers
  dat$q1_a <- case_when(dat$q1_a == "Fund A" ~ "Mid Risk",
                        dat$q1_a == "Fund B" ~ "More Risk",
                        dat$q1_a == "Fund C" ~ "Less Risk")
  
  dat$q2_a <- case_when(dat$q2_a == "Fund A" ~ "Mid Risk",
                        dat$q2_a == "Fund B" ~ "More Risk",
                        dat$q2_a == "Fund C" ~ "Less Risk")
  
  dat$q3_a <- case_when(dat$q3_a == "Fund A" ~ "Mid Risk",
                        dat$q3_a == "Fund B" ~ "More Risk",
                        dat$q3_a == "Fund C" ~ "Less Risk")
  
  dat$q4_a <- case_when(dat$q4_a == "Fund A" ~ "TDF",
                        dat$q4_a == "Fund B" ~ "Floor",
                        dat$q4_a == "Fund C" ~ "Buffer")
  
  
  dat$group_assigned <- case_when(dat$group_assigned == 1 ~ "1. Base",
                                  dat$group_assigned == 2 ~ "2. Default Option",
                                  dat$group_assigned == 3 ~ "3. Distributional Info",
                                  dat$group_assigned == 4 ~ "4. Dist. and Default")
  
  (log(.8)-log(.4))*.12
  
  
  
  btab <- dat %>% 
    mutate(Male = as.numeric(aq1 == "Male"),
           White = as.numeric(aq2 == "Caucasian – not Hispanic" |
                                aq2 == "Caucasian â€“ not Hispanic"),
           Married = as.numeric(aq3 == "Married"),
           HS = as.numeric(aq4 == "High school graduate or GED equivalent")) %>% 
    group_by(group_assigned) %>% 
    summarize(n = n(),
              age_current_mean = mean(age_current),
              age_current_sd = sd(age_current),
              salary_mean = mean(salary_current/1000),
              salary_sd = sd(salary_current/1000),
              Male_mean = mean(Male)*100,
              Male_sd = sd(Male)*100,
              White_mean = mean(White)*100,
              White_sd = sd(White)*100,
              Married_mean = mean(Married)*100,
              Married_sd = sd(Married)*100,
              HS_mean = mean(HS)*100,
              HS_sd = sd(HS)*100)
  
  
  
  # Balance Table
  mdif <- function(m1, m2, s1, s2, n1, n2)
  {
    (m1 - m2)/(sqrt(s1/n1 + s2/n2))
  }
  
  # Formatting
  f <- function(x, digits = 2){
    if(x>100){  format(x, digits = digits, nsmall = 0, big.mark = ",") 
    } else {
      format(x, digits = digits, nsmall = 2, big.mark = ",") 
    }
  }
  
  cat("Current Age & ",
      f(btab$age_current_mean[1]), "&",
      f(btab$age_current_mean[2]), "&",
      f(btab$age_current_mean[3]), "&",
      f(btab$age_current_mean[4]), "\\\\ & (",
      f(btab$age_current_sd[1]), ") & (",
      f(btab$age_current_sd[2]), ") & (",
      f(btab$age_current_sd[3]), ") & (",
      f(btab$age_current_sd[4]), ") \\\\ \\\\
      Average Salary (\\$000s) & ",
      f(btab$salary_mean[1]), "&",
      f(btab$salary_mean[2]), "&",
      f(btab$salary_mean[3]), "&",
      f(btab$salary_mean[4]), "\\\\ & (",
      f(btab$salary_sd[1]), ") & (",
      f(btab$salary_sd[2]), ") & (",
      f(btab$salary_sd[3]), ") & (",
      f(btab$salary_sd[4]), ") \\\\ \\\\
      Percent Male & ",
      f(btab$Male_mean[1]), "&",
      f(btab$Male_mean[2]), "&",
      f(btab$Male_mean[3]), "&",
      f(btab$Male_mean[4]), "\\\\ & (",
      f(btab$Male_sd[1]), ") & (",
      f(btab$Male_sd[2]), ") & (",
      f(btab$Male_sd[3]), ") & (",
      f(btab$Male_sd[4]), ") \\\\ \\\\
      Percent White & ",
      f(btab$White_mean[1]), "&",
      f(btab$White_mean[2]), "&",
      f(btab$White_mean[3]), "&",
      f(btab$White_mean[4]), "\\\\ & (",
      f(btab$White_sd[1]), ") & (",
      f(btab$White_sd[2]), ") & (",
      f(btab$White_sd[3]), ") & (",
      f(btab$White_sd[4]), ") \\\\ \\\\
      Percent Married & ",
      f(btab$Married_mean[1]), "&",
      f(btab$Married_mean[2]), "&",
      f(btab$Married_mean[3]), "&",
      f(btab$Married_mean[4]), "\\\\& (",
      f(btab$Married_sd[1]), ") & (" ,
      f(btab$Married_sd[2]), ") & (",
      f(btab$Married_sd[3]), ") & (",
      f(btab$Married_sd[4]), ") \\\\ \\\\
      Percent No College & ",
      f(btab$HS_mean[1]), "&",
      f(btab$HS_mean[2]), "&",
      f(btab$HS_mean[3]), "&",
      f(btab$HS_mean[4]), "\\\\& (",
      f(btab$HS_sd[1]), ") & (",
      f(btab$HS_sd[2]), ") & (",
      f(btab$HS_sd[3]), ") & (",
      f(btab$HS_sd[4]), ") \\\\
      \\bottomrule \\\\
      Number of Observations & \\multicolumn{1}{c}{",
      f(btab$n[1]), "} & \\multicolumn{1}{c}{",
      f(btab$n[2]), "} & \\multicolumn{1}{c}{",
      f(btab$n[3]), "} & \\multicolumn{1}{c}{",
      f(btab$n[4]), "}\\\\", sep="",
      file = paste0(opath, "Tables/bal_table_existing_RA.tex"))
  
  
  
  
  
  # Financial Literacy ------------------------------------------------------
  dat$fin_lit <- if_else(dat$aq8 %in% c("Somewhat sure â€“ I understand most of what I will need to know", "Very sure â€“ I understand money management very well"),
                         1, 0)
  
  
  # subset dat to just group 1 and 2
  dat12 <- 
    dat %>% 
    filter(group_assigned %in% c(1,2))
  
  
  
  
  
  
  
  
  # Data Manipulation -------------------------------------------------------
  dat$risk_aversion <- ifelse(grepl("200", dat$la7), "Less Risk Averse", "More Risk Averse")
  summary(dat$risk_aversion)
  
  dat$loss_aversion <- 
    ifelse(
      (
        (
          grepl("sure", dat$aq10) &
            grepl("probability", dat$aq11) 
        ) |
          (
            grepl("sure", dat$aq12) &
              grepl("probability", dat$aq13)
          )
      ), 
      "More Loss Averse",
      "Less Loss Averse")
  
  table(dat$loss_aversion)
  
  dat$time <- dat$finish - dat$start
  
  dat$marg_rate_1 <- ifelse(is.na(dat$q1_c1), dat$q1_c2, dat$q1_c1)
  summary(dat$marg_rate_1)
  
  dat$marg_rate_2 <- ifelse(is.na(dat$q2_c1), dat$q2_c2, dat$q2_c1)
  summary(dat$marg_rate_2)
  
  dat$marg_rate_3 <- ifelse(is.na(dat$q3_c1), dat$q3_c2, dat$q3_c1)
  summary(dat$marg_rate_3)
  
  dat$marg_rate_4 <- ifelse(is.na(dat$q4_c1), dat$q4_c2, dat$q4_c1)
  summary(dat$marg_rate_4)
  
  # Binary Versions
  dat$retire_time_bin <- ifelse((dat$age_retire - dat$age_current) > median(dat$age_retire - dat$age_current),
                                "Near Retirement",
                                "Not Near Retirement")
  dat$retire_size_bin <- ifelse((dat$retire_plan) > median(dat$retire_plan),
                                "Large Current Plan",
                                "Small Current Plan")
  
  
  
  summary(dat$change1)
  summary(dat$change2)
  summary(dat$change3)
  summary(dat$change4)
  
  .16*(log(.4)- log(.2))/.3
  
  
  
  dat$gross_expense_ratio_random <- case_when(dat$gross_expense_ratio_random == 2 ~ .003,
                                              dat$gross_expense_ratio_random == 3 ~ .005,
                                              dat$gross_expense_ratio_random == 4 ~ .007)
  
  a <- lm(I(1- dat$change1) ~ dat$group_assigned:I(log(dat$gross_expense_ratio_random))) %>% summary()
  b <- lm(I(1- dat$change2) ~ dat$group_assigned:I(log(dat$gross_expense_ratio_random))) %>% summary()
  c <- lm(I(1- dat$change3) ~ dat$group_assigned:I(log(dat$gross_expense_ratio_random))) %>% summary()
  d <- lm(I(1- dat$change4) ~ dat$group_assigned:I(log(dat$gross_expense_ratio_random))) %>% summary()
  
  cat("log(Expense Ratio):Base Group & ", 
      f(a$coefficients[2,1], 3), " & ",
      f(b$coefficients[2,1], 3), " & ",
      f(c$coefficients[2,1], 3), " & ",
      f(d$coefficients[2,1], 3), " \\\\ 
  & (", 
      f(a$coefficients[2,3], 3), ") & (",
      f(b$coefficients[2,3], 3), ") & (",
      f(c$coefficients[2,3], 3), ") & (",
      f(d$coefficients[2,3], 3), ") \\\\ \\\\
  log(Expense Ratio):Default Group& ", 
      f(a$coefficients[3,1], 3), " & ",
      f(b$coefficients[3,1], 3), " & ",
      f(c$coefficients[3,1], 3), " & ",
      f(d$coefficients[3,1], 3), " \\\\  
  & (", 
      f(a$coefficients[3,3], 3), ") & (",
      f(b$coefficients[3,3], 3), ") & (",
      f(c$coefficients[3,3], 3), ") & (",
      f(d$coefficients[3,3], 3), ") \\\\ \\\\
  log(Expense Ratio):Distribution Group & ", 
      f(a$coefficients[4,1], 3), " & ",
      f(b$coefficients[4,1], 3), " & ",
      f(c$coefficients[4,1], 3), " & ",
      f(d$coefficients[4,1], 3), " \\\\ 
  & (",  
      f(a$coefficients[4,3], 3), ") & (",
      f(b$coefficients[4,3], 3), ") & (",
      f(c$coefficients[4,3], 3), ") & (",
      f(d$coefficients[4,3], 3), ") \\\\ \\\\
  log(Expense Ratio):Default and Distribution & ", 
      f(a$coefficients[5,1], 3), " & ",
      f(b$coefficients[5,1], 3), " & ",
      f(c$coefficients[5,1], 3), " & ",
      f(d$coefficients[5,1], 3), " \\\\ 
  & (",
      f(a$coefficients[5,3], 3), ") & (",
      f(b$coefficients[5,3], 3), ") & (",
      f(c$coefficients[5,3], 3), ") & (",
      f(d$coefficients[5,3], 3), ") \\\\", file = paste0(opath, "tables/switching_existing_RA.tex"))
  
      
      
      
      
      
      
      
      
      
      
      
      
      
  
  l <- which(grepl("Dist", dat$group_assigned))
  
  lm(change1 ~ I(log(gross_expense_ratio_random))*loss_aversion + I(log(gross_expense_ratio_random))*risk_aversion, data = dat) %>% summary()
  lm(change2 ~ I(log(gross_expense_ratio_random))*loss_aversion + I(log(gross_expense_ratio_random))*risk_aversion, data = dat) %>% summary()
  lm(change3 ~ I(log(gross_expense_ratio_random))*loss_aversion + I(log(gross_expense_ratio_random))*risk_aversion, data = dat) %>% summary()
  lm(change4 ~ I(log(gross_expense_ratio_random))*loss_aversion + I(log(gross_expense_ratio_random))*risk_aversion, data = dat) %>% summary()
  
  
  lm(I(q4_a == "Buffer RILA")  ~ group_assigned*loss_aversion, data = dat) %>% summary()
  
  # Initial Choice ----------------------------------------------------------
  # Reshape the data
  dat_long <- dat %>%
    # Keep only the columns we're interested in for the plot
    select(q1_a, q2_a, q3_a) %>%
    # Convert to long format
    gather(key = "question_type", value = "value")
  
  # Create the plot
  p <- dat_long %>%
    ggplot(aes(x = value)) +
    geom_bar() +
    # Map facets to the question type
    facet_wrap(~ question_type, scales = "free", labeller = as_labeller(c(q1_a = "TDF - Initial Choice", q2_a = "TD-RILA with Floor - Initial Choice", q3_a = "TD-RILA with Buffer - Initial Choice"))) +
    ylab("Count") +
    xlab("Investment Choice - Risk Level") +
    theme(
      text = element_text(size = t1),  # Global text size
      axis.title = element_text(size = t2),  # Axis title size
      axis.text = element_text(size = t2)  # Axis text size
    )
  
  # Display or save the plot
  if(sa == 1) {
    ggsave(file = paste0(opath, "Figures/combined_plot_existing_RA.jpg"), plot = p, width = width, height = height)
  } else {
    print(p)
  }
  
  
  
  # q1a
  dat %>% 
    # group_by(group_assigned) %>% 
    # mutate(count = n()) %>% 
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q1_a)) +
    geom_bar() +
    # facet_grid(. ~ group_assigned) + 
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Count") + 
    xlab("TDF - Initial Choice") + 
    theme(
      text = element_text(size = t1), # Global text size
      axis.title = element_text(size =t2), # Axis title size
      axis.text = element_text(size =t2) # Axis text size
    )
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q1a1_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  dat %>% 
    # group_by(group_assigned) %>% 
    # mutate(count = n()) %>% 
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q2_a)) +
    geom_bar() +
    # facet_grid(. ~ group_assigned) + 
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Count") + 
    xlab("TD-RILA with Floor - Initial Choice") + 
    theme(
      text = element_text(size = t1), # Global text size
      axis.title = element_text(size =t2), # Axis title size
      axis.text = element_text(size =t2) # Axis text size
    )
  
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q2a1_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  
  dat %>% 
    # group_by(group_assigned) %>% 
    # mutate(count = n()) %>% 
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q3_a)) +
    geom_bar() +
    # facet_grid(. ~ group_assigned) + 
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Count") + 
    xlab("TD-RILA with Buffer - Initial Choice") + 
    theme(
      text = element_text(size = t1), # Global text size
      axis.title = element_text(size =t2), # Axis title size
      axis.text = element_text(size =t2) # Axis text size
    )
  
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q3a1_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  dat %>% 
    filter(group_assigned == "1. Base") %>%
    # mutate(count = n()) %>% 
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q4_a)) +
    geom_bar() +
    # facet_grid(. ~ group_assigned) + 
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Count") + 
    xlab("TDF - Initial Choice") + 
    theme(
      text = element_text(size = t1), # Global text size
      axis.title = element_text(size =t2), # Axis title size
      axis.text = element_text(size =t2) # Axis text size
    )
  
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q4a1_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  
  # q1a
  dat %>% 
    group_by(group_assigned) %>% 
    mutate(count = n()) %>% 
    group_by(group_assigned, q1_a) %>%
    summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q1_a, y = perc)) +
    geom_col() +
    facet_grid(. ~ group_assigned) + 
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") + 
    xlab("TDF - Initial Choice") + 
    theme(
      text = element_text(size = t1), # Global text size
      axis.title = element_text(size =t2), # Axis title size
      axis.text = element_text(size =t2) # Axis text size
    )
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q1a_existing_RA.jpg"),
                     width = width2,
                     height = height)
  
  # q2a
  dat %>% 
    group_by(group_assigned) %>% 
    mutate(count = n()) %>% 
    group_by(group_assigned, q2_a) %>%
    summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q2_a, y = perc)) +
    geom_col() +
    facet_grid(. ~ group_assigned)+ 
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") + 
    xlab("Floor RILA - Initial Choice") + 
    theme(
      text = element_text(size = t1), # Global text size
      axis.title = element_text(size =t2), # Axis title size
      axis.text = element_text(size =t2) # Axis text size
    )
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q2a_existing_RA.jpg"),
                     width = width2,
                     height = height)
  
  # q3a
  dat %>% 
    group_by(group_assigned) %>% 
    mutate(count = n()) %>% 
    group_by(group_assigned, q3_a) %>%
    summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q3_a, y = perc)) +
    geom_col() +
    facet_grid(. ~ group_assigned)+ 
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") + 
    xlab("Buffer RILA - Initial Choice") + 
    theme(
      text = element_text(size = t1), # Global text size
      axis.title = element_text(size =t2), # Axis title size
      axis.text = element_text(size =t2) # Axis text size
    )
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q3a_existing_RA.jpg"),
                     width = width2,
                     height = height)
  
  # q4a
  dat %>% 
    group_by(group_assigned) %>% 
    mutate(count = n()) %>% 
    group_by(group_assigned, q4_a) %>%
    summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(. ~ group_assigned)+ 
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") + 
    xlab("Combo - Initial Choice")
  
  # If you want to save run this first :
  # sa <- 1
  
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q4a_existing_RA.jpg"),
                     width = width2,
                     height = height)
  
  # q4a - ra
  dat %>% 
    group_by(group_assigned, risk_aversion) %>% 
    mutate(count = n()) %>% 
    group_by(group_assigned, risk_aversion, q4_a) %>%
    summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(risk_aversion ~ group_assigned) + 
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") + 
    xlab("Combo - Initial Choice")
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q4a_ra_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  
  # q4a - la
  dat %>% 
    group_by(group_assigned, loss_aversion) %>% 
    mutate(count = n()) %>% 
    group_by(group_assigned, loss_aversion, q4_a) %>%
    summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(loss_aversion ~ group_assigned) + 
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") + 
    xlab("Combo - Initial Choice")
  
  # q4a - la
  dat %>% 
    group_by(group_assigned, loss_aversion) %>% 
    mutate(count = n()) %>% 
    group_by(group_assigned, loss_aversion, q4_a) %>%
    summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(loss_aversion ~ group_assigned) + 
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") + 
    xlab("Combo - Initial Choice")
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q4a_la_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  # q4a - la
  dat %>% 
    group_by(group_assigned, retire_time_bin) %>% 
    mutate(count = n()) %>% 
    group_by(group_assigned, retire_time_bin, q4_a) %>%
    summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(retire_time_bin ~ group_assigned) + 
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") + 
    xlab("Combo - Initial Choice")
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q4a_rt_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  
  # q4a - la
  dat %>% 
    group_by(group_assigned, retire_size_bin) %>% 
    mutate(count = n()) %>% 
    group_by(group_assigned, retire_size_bin, q4_a) %>%
    summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(retire_size_bin ~ group_assigned) + 
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") + 
    xlab("Combo - Initial Choice")
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q4a_rs_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  
    dat$fin_lit <- if_else(dat$aq16 %in%  c("More than today", "Exactly the same"), 0, 1)
    
  dat %>% 
    group_by(group_assigned, fin_lit) %>% 
    mutate(count = n()) %>% 
    group_by(group_assigned, fin_lit, q4_a) %>%
    summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(fin_lit ~ group_assigned) + 
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") + 
    xlab("Combo - Initial Choice")
  
  # Demand Elasticities ------------------------------------------
  
  # q1a
  dat %>% 
    group_by(group_assigned) %>% 
    # mutate(count = n()) %>% 
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = marg_rate_1)) +
    geom_density() +
    facet_grid(.~ group_assigned) + 
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Density") + 
    xlab("TDF - Marginal Fee")
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q1_marg_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  # q1b
  dat %>% 
    group_by(group_assigned) %>% 
    # mutate(count = n()) %>% 
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = marg_rate_2)) +
    geom_density() +
    facet_grid(. ~ group_assigned) + 
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Density") + 
    xlab("Floor RILA - Marginal Fee")
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q2_marg_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  # q1a
  dat %>% 
    group_by(group_assigned) %>% 
    # mutate(count = n()) %>% 
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = marg_rate_3)) +
    geom_density() +
    facet_grid(. ~ group_assigned) + 
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Density") + 
    xlab("Buffer RILA - Marginal Fee")
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q3_marg_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  # q4a
  dat %>% 
    group_by(group_assigned) %>% 
    # mutate(count = n()) %>% 
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = marg_rate_4)) +
    geom_density() +
    facet_grid(.  ~ group_assigned) + 
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Density") + 
    xlab("Crossplan - Marginal Fee")
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q4_marg_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  # q4a
  dat %>% 
    group_by(group_assigned, risk_aversion) %>% 
    # mutate(count = n()) %>% 
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = marg_rate_4)) +
    geom_density() +
    facet_grid(risk_aversion  ~ group_assigned) + 
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Density") + 
    xlab("Crossplan - Marginal Fee")
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q4_marg_ra_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  # q4a
  dat %>% 
    group_by(group_assigned, loss_aversion) %>% 
    # mutate(count = n()) %>% 
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = marg_rate_4)) +
    geom_density() +
    facet_grid(loss_aversion  ~ group_assigned) + 
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Density") + 
    xlab("Crossplan - Marginal Fee")
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q4_marg_la_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  # q4a
  dat %>% 
    group_by(group_assigned, retire_time_bin) %>% 
    # mutate(count = n()) %>% 
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = marg_rate_4)) +
    geom_density() +
    facet_grid(retire_time_bin  ~ group_assigned) + 
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Density") + 
    xlab("Crossplan - Marginal Fee")
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q4_marg_rt_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  # q4a
  dat %>% 
    group_by(group_assigned, fin_lit) %>% 
    # mutate(count = n()) %>% 
    # group_by(group_assigned, q1_a) %>%
    # summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = marg_rate_4)) +
    geom_density() +
    facet_grid(fin_lit  ~ group_assigned) + 
    # scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Density") + 
    xlab("Crossplan - Marginal Fee")
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q4_marg_rs_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  # Income-based analysis
  dat$income_group <- ifelse(dat$salary_current > median(dat$salary_current), 
                             "High Income", "Low Income")
  
  # Plot for income groups
  dat %>% 
    group_by(group_assigned, income_group) %>% 
    mutate(count = n()) %>% 
    group_by(group_assigned, income_group, q4_a) %>%
    summarise(perc = n()/count) %>%
    distinct() %>% 
    ggplot(aes(x = q4_a, y = perc)) +
    geom_col() +
    facet_grid(income_group ~ group_assigned) + 
    scale_y_continuous(labels = scales::percent_format(scale = 100)) +
    ylab("Percent") + 
    xlab("Combo - Initial Choice") +
    theme(
      text = element_text(size = t1), # Global text size
      axis.title = element_text(size = t2), # Axis title size
      axis.text = element_text(size = t2) # Axis text size
    )
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q4a_income_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  # Marginal rate analysis for income groups
  dat %>% 
    group_by(group_assigned, income_group) %>% 
    distinct() %>% 
    ggplot(aes(x = marg_rate_4)) +
    geom_density() +
    facet_grid(income_group ~ group_assigned) + 
    ylab("Density") + 
    xlab("Crossplan - Marginal Fee") +
    theme(
      text = element_text(size = t1), # Global text size
      axis.title = element_text(size = t2), # Axis title size
      axis.text = element_text(size = t2) # Axis text size
    )
  
  if(sa == 1) ggsave(file = paste0(opath, "Figures/q4_marg_income_existing_RA.jpg"),
                     width = width,
                     height = height)
  
  # Statistical test for differences in initial choices
  income_choice_test <- chisq.test(table(dat$income_group, dat$q4_a))
  
  # Statistical test for differences in marginal rates
  income_marg_test <- t.test(marg_rate_4 ~ income_group, data = dat)
  
  # Print test results
  cat("Chi-square test for initial choices by income group:\n")
  print(income_choice_test)
  
  cat("\nt-test for marginal rates by income group:\n")
  print(income_marg_test)
  
  # Function to calculate confidence interval
  calc_ci <- function(x, n) {
    p <- mean(x)
    se <- sqrt(p * (1 - p) / n)
    ci <- 1.96 * se
    return(c(p - ci, p + ci))
  }
  
  # Function to prepare data for plotting
  prepare_plot_data <- function(data, question) {
    data %>%
      group_by(group_assigned) %>%
      summarize(
        Less_Risk = mean(!!sym(question) == "Less Risk"),
        Mid_Risk = mean(!!sym(question) == "Mid Risk"),
        More_Risk = mean(!!sym(question) == "More Risk"),
        n = n(),
        Mid_Risk_CI = list(calc_ci(!!sym(question) == "Mid Risk", n())),
        More_Risk_CI = list(calc_ci(!!sym(question) == "More Risk", n()))
      ) %>%
      mutate(
        Mid_Risk_Diff = Mid_Risk - Less_Risk,
        More_Risk_Diff = More_Risk - Less_Risk,
        Mid_Risk_CI_Low = Mid_Risk_Diff - (Mid_Risk_CI[[1]][2] - Mid_Risk_CI[[1]][1]) / 2,
        Mid_Risk_CI_High = Mid_Risk_Diff + (Mid_Risk_CI[[1]][2] - Mid_Risk_CI[[1]][1]) / 2,
        More_Risk_CI_Low = More_Risk_Diff - (More_Risk_CI[[1]][2] - More_Risk_CI[[1]][1]) / 2,
        More_Risk_CI_High = More_Risk_Diff + (More_Risk_CI[[1]][2] - More_Risk_CI[[1]][1]) / 2
      )
  }
  
  # Prepare data for each question
  q1_data <- prepare_plot_data(dat, "q1_a")
  q2_data <- prepare_plot_data(dat, "q2_a")
  q3_data <- prepare_plot_data(dat, "q3_a")
  
  # Function to create plot
  create_plot <- function(data, title) {
    ggplot(data, aes(x = group_assigned)) +
      geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
      geom_pointrange(aes(y = Mid_Risk_Diff, ymin = Mid_Risk_CI_Low, ymax = Mid_Risk_CI_High, color = "Mid Risk"), 
                     position = position_dodge(width = 0.5)) +
      geom_pointrange(aes(y = More_Risk_Diff, ymin = More_Risk_CI_Low, ymax = More_Risk_CI_High, color = "More Risk"), 
                     position = position_dodge(width = 0.5)) +
      scale_color_manual(values = c("Mid Risk" = "blue", "More Risk" = "red")) +
      labs(title = title,
           x = "Treatment Group",
           y = "Difference in Proportion (Relative to Less Risk)",
           color = "Risk Level") +
      theme_minimal() +
      theme(
        text = element_text(size = t1),
        axis.title = element_text(size = t2),
        axis.text = element_text(size = t2),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )
  }
  
  # Create plots
  p1 <- create_plot(q1_data, "TDF - Initial Choice")
  p2 <- create_plot(q2_data, "TD-RILA with Floor - Initial Choice")
  p3 <- create_plot(q3_data, "TD-RILA with Buffer - Initial Choice")
  
  # Save plots
  if(sa == 1) {
    ggsave(file = paste0(opath, "Figures/q1a_relative_existing_RA.jpg"), plot = p1, width = width2, height = height)
    ggsave(file = paste0(opath, "Figures/q2a_relative_existing_RA.jpg"), plot = p2, width = width2, height = height)
    ggsave(file = paste0(opath, "Figures/q3a_relative_existing_RA.jpg"), plot = p3, width = width2, height = height)
  }
  
  # For q4_a (Combo), we need to handle it differently as it has different categories
  q4_data <- dat %>%
    group_by(group_assigned) %>%
    summarize(
      TDF = mean(q4_a == "TDF"),
      Floor = mean(q4_a == "Floor"),
      Buffer = mean(q4_a == "Buffer"),
      n = n(),
      Floor_CI = list(calc_ci(q4_a == "Floor", n())),
      Buffer_CI = list(calc_ci(q4_a == "Buffer", n()))
    ) %>%
    mutate(
      Floor_Diff = Floor - TDF,
      Buffer_Diff = Buffer - TDF,
      Floor_CI_Low = Floor_Diff - (Floor_CI[[1]][2] - Floor_CI[[1]][1]) / 2,
      Floor_CI_High = Floor_Diff + (Floor_CI[[1]][2] - Floor_CI[[1]][1]) / 2,
      Buffer_CI_Low = Buffer_Diff - (Buffer_CI[[1]][2] - Buffer_CI[[1]][1]) / 2,
      Buffer_CI_High = Buffer_Diff + (Buffer_CI[[1]][2] - Buffer_CI[[1]][1]) / 2
    )
  
  p4 <- ggplot(q4_data, aes(x = group_assigned)) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
    geom_pointrange(aes(y = Floor_Diff, ymin = Floor_CI_Low, ymax = Floor_CI_High, color = "Floor RILA"), 
                    position = position_dodge(width = 0.5)) +
    geom_pointrange(aes(y = Buffer_Diff, ymin = Buffer_CI_Low, ymax = Buffer_CI_High, color = "Buffer RILA"), 
                    position = position_dodge(width = 0.5)) +
    scale_color_manual(values = c("Floor RILA" = "blue", "Buffer RILA" = "red")) +
    labs(title = "Combo - Initial Choice",
         x = "Treatment Group",
         y = "Difference in Proportion (Relative to TDF)",
         color = "RILA Type") +
    theme_minimal() +
    theme(
      text = element_text(size = t1),
      axis.title = element_text(size = t2),
      axis.text = element_text(size = t2),
      axis.text.x = element_text(angle = 45, hjust = 1)
    )
  
  if(sa == 1) {
    ggsave(file = paste0(opath, "Figures/q4a_relative_existing_RA.jpg"), plot = p4, width = width2, height = height)
  }